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ABSTRACT 

A recipe is presented to construct an analytic, self-consistent model of a non-barotropic 
neutron star with a poloidal-toroidal field of arbitrary multipole order, whose toroidal 
component is confined in a torus around the neutral curve inside the star, as in numer¬ 
ical simulations of twisted tori. The recipe takes advantage of magnetic-field-aligned 
coordinates to ensure continuity of the mass density at the surface of the torus. The 
density perturbation and ellipticity of such a star are calculated in general and for 
the special case of a mixed dipole-quadrupole field as a worked example. The calcula¬ 
tion generalises previous work restricted to dipolar, poloidal-toroidal and multipolar, 
poloidal-only configurations. The results are applied, as an example, to magnetars 
whose observations (e.g., spectral features and pulse modulation) indicate that the 
internal magnetic fields may be at least one order of magnitude stronger than the 
external fields, as inferred from their spin downs, and are not purely dipolar. 

Key -words: MHD - stars: magnetic field - stars: interiors - stars: neutron - gravi¬ 
tational waves 


1 INTRODUCTION 

Neutron stars, particularly the subset known as magnetars, possess some of the strongest known, naturally occurring magnetic 
fields. The external field strength of a neutron star is inferred from its spin down. The internal field cannot be measured 
directly. Observations of bursts and giant flares (loka 2001; Corsi & Owen 2011) and precession (Makishima et al. 2014) in 
magnetars have been interpreted to indicate that the internal field exceeds the external field by at least one order of magnitude. 
Numerical simulations favour a ‘twisted torus’ magnetic configuration (Braithwaite & Nordlund 2006; Braithwaite & Spruit 
2006). 

Magnetic stresses deform a neutron star (Chandrasekhar & Fermi 1953; Ferraro 1954; Goosens 1972; Katz 1989; Cutler 
2002; Payne & Melatos 2004; Haskell et al. 2008; Mastrano et al. 2011), leading to detectable levels of gravitational radiation 
under certain conditions (Bonazzola & Gourgoulhon 1996; Melatos & Payne 2005; Stella et al. 2005; Haskell et al. 2008; 
DairOsso, Shore, & Stella 2009). Upper limits from gravitational wave non-detections to date can be used to set upper limits 
on the star’s ellipticity, and hence to constrain the strength and topology of the internal field in principle (Gutler 2002; 
DaU’Osso, Shore, & Stella 2009; Mastrano et al. 2011; Mastrano & Melatos 2012). 

Neutron star magnetic fields are approximately dipolar at radio emission altitudes (Lyne & Manchester 1988; Chung & 
Melatos 2011; Burnett & Melatos 2014) and in the outer magnetosphere, where high-energy emissions originate (Romani & 
Yadigaroglu 1995; Lyutikov, Otte, & McCann 2012). However, some observations can be interpreted as signatures of higher- 
order multipoles close to the surface. Examples include intermittent radio emission from pulsars beyond the pair-cascade 
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‘death line’ (Young, Manchester, & Johnston 1999; Camilo et al. 2000; Gil & Mitra 2001; McLaughlin et al. 2003; Medin & 
Lai 2010), the pulse profile of SGR 1900+14 following its 1998 August 27 giant flare (Feroci et al. 2001; Thompson & Duncan 
2001; Thompson, Lyutikov, & Kulkarni 2002), cyclotron resonant scattering line energies of some accretion-powered X-ray 
pulsars (Nishimura 2005; Priymak, Melatos, & Lasky 2014), the anomalous braking index of some radio pulsars (Barsukov & 
Tsygan 2010), and X-ray spectral features of SGR 0418+5729 (Giiver, Ozel F., & G6gu§ 2008; Gfiver, G6gu§, & Ozel 2011; 
Tiengo et al. 2013) and PSR J0821—4300 (Gotthelf, Halpern, & Alford 2013) . While the dipole component of the magnetic 
field is readily inferred from its spin down, the putative higher-order multipoles do not contribute much to the spin-down rate. 

In analytically studying higher-order multipoles through (say) their gravitational wave emission, one confronts a key 
technical challenge, namely ensuring that the magnetically-induced pressure and density perturbations are continuous every¬ 
where. If the field has both poloidal and toroidal components, and the form of the toroidal component is not chosen wisely, one 
ends up with discontinuous pressure and density proHles. The ‘correct’ form of the toroidal held (i.e., one that is analytically 
tractable and does not lead to discontinuities) is not immediately obvious. In previous papers, we investigated the relation 
between held conhguration and stellar ellipticity in a stratihed, non-barotropic star with an axisymmetric, purely dipolar, 
poloidal and toroidal held (Mastrano et al. 2011) and an axisymmetric, multipolar, purely poloidal held (Mastrano, Lasky, 
& Melatos 2013). Realistic neutron stars are non-barotropic; they are stably stratihed by a composition gradient (Pethick 
1992; Reisenegger & Goldreich 1992; Reisenegger 2009). Because a non-barotropic star allows arbitrary relative poloidal and 
toroidal held strengths, it can easily accommodate the strong internal toroidal helds suggested by observations (loka 2001; 
Corsi & Owen 2011; Makishima et al. 2014) and simulations (Braithwaite & Nordlund 2006). However, like other authors, 
we were unable to construct a physically and mathematically consistent poloidal-toroidal multipolar held. We found that the 
simple forms of an axisymmetric held used by Mastrano et al. (2011) and others lead to unphysical density discontinuities 
at the boundary of the toroidal held region. There appears to be no record in the literature of an analytic poloidal-toroidal 
construction for multipoles of higher order than dipole. This is unfortunate, because once this issue is dealt with, linear 
combinations of multipoles can generate realistic analytic models which complement numerical simulations helpfully. 

In this short methods paper, we present a new method for constructing such a consistent axisymmetric poloidal-toroidal 
held analytically for any multipole. In Sec. 2, we recap the results of Mastrano et al. (2011) and Mastrano, Lasky, & Melatos 
(2013), and we describe the method for calculating e for a generalised multipolar held. In the Appendix, we show why the new 
method presented here is necessary, by way of an explicit example, where the old methods fail. In Sec. 3, we present a worked 
example of a mixed dipole-quadrupole, poloidal-toroidal held, constructed with the new method. We give explicit expressions 
for the held and calculate the deformation due to this held. In Sec. 4, we apply the results briehy to the magnetars SGR 
0418+5729 and 4U 0142+61 to give a havour of the sort of astrophysical questions that the method can help to study. Lastly, 
in Sec. 5, we summarize our hndings. We emphasize that our analytic held models should complement, rather than supplant, 
large-scale numerical simulations. They are not meant to be true representations of neutron star helds but rather to provide 
helpful toy models, which are easy to manipulate and whose general properties are easy to parametrize and understand. 


2 MASS ELLIPTICITY 

In this section, we describe a general method for constructing an analytic model of a non-barotropic star in hydromagnetic 
equilibrium, where the poloidal component of the magnetic held is a linear combination of multipolar helds, and the toroidal 
component is located around the neutral curve (the circumstellar curve inside the star where the poloidal component vanishes). 
This method is needed to ensure continuity of the density held at the boundary of the magnetic torus; an example of what 
can go wrong without it is given in the Appendix. The new feature of the method is a held-aligned coordinate transformation 
which allows one to solve the force balance equation exactly in closed form for the density perturbations, generated by a 
poloidal magnetic stream function of arbitrary multipole order. Previous calculations (Mastrano et al. 2011; Mastrano & 
Melatos 2012; Yoshida 2013; Dall’Osso et al. 2014) are restricted to dipolar poloidal helds. 


2.1 Hydromagnetic force balance 

We decompose the magnetic held into its poloidal and toroidal components and express it in dimensionless spherical polar 
coordinates (r, 6, <j)) in the usual way (Ghandrasekhar 1956; Mastrano et al. 2011; Mastrano & Melatos 2012; Mastrano, Lasky, 
& Melatos 2013), viz. 


B = Bolrip\7a{r,e) X \7(f) + ritl3{a)\7(f)], (1) 

where Bq parametrizes the overall strength of the held, r/p and rjt set the relative strengths of the poloidal and toroidal 
components respectively {rjp = 1 without loss of generality), a{r,6) is the poloidal magnetic stream function, and the function 
/3(a) dehnes the toroidal held component. In general, the stream function a{r,0) can be expanded as a linear combination of 
multipolar stream functions ai{r,0) up to any desired order n, viz. 
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n 

a{r,e) = '^aiai{r,e), ( 2 ) 

i=i 

where the subscripts denote the multipole orders and ai is a dimensionless weight. We consider separable stream functions of 
the form ai{r,6) = fi{r)gi{9) in this paper. The function p must be a function of a to ensure that the magnetic force has no 
azimuthal component, which cannot be balanced in magnetohydrostatic equilibrium given a field of the form (1) (Mastrano 
et al. 2011; Mastrano, Lasky, & Melatos 2013). 

The magnetic energy density is < 10“® of the gravitational energy density, even in magnetars. Therefore, we can treat 
the magnetic force as a perturbation on a background hydrostatic equilibrium and write the hydromagnetic force balance 
equation as 


— (V X B) X B = Vap + ^pV-h, (3) 

po 

to first order in /{fiop) in the Cowling approximation (5>1> = 0), where p is the zeroth-order pressure, p is the zeroth-order 
density, <1? is the gravitational potential, and Sp, Sp, J"!? are perturbations of the latter three quantities. Because we do not 
assume a barotropic star, the density perturbation Sp does not have to be a function solely of the pressure perturbation Sp, 
and therefore the equation of state imposes no restrictions on the field structure. Physically, this means that the imposed 
magnetic field sets the density and pressure perturbations, but the resulting perturbations do not restrict the magnetic field 
in turn. Therefore, unlike some previous works [e.g., Haskell et al. (2008); Lander & Jones (2009); Ciolfi, Ferrari, & Gualtieri 
(2010)], we do not specify a barotropic equation of state and then solve the Grad-Shafranov equation for the magnetic field 
configuration. Instead, we specify the magnetic field whose effects we wish to investigate [to be determined from observations, 
e.g., neutron star spin down, gravitational wave upper limits, and radio polarization (Chung & Melatos 2011; Burnett & 
Melatos 2014)] then calculate the density perturbations that the field causes. 

We characterize the magnetic deformation of the star by its ellipticity e. 




lo ’ 

where 7o is the moment of inertia of the unperturbed spherical star, the moment-of-inertia tensor is given by 


( 4 ) 


Ijk ^ RI [ A^x[p{r) + Sp{r,e)]{r^Sjk - XjXk), (5) 

Jv 

R» is the stellar radius, and the integral is taken over the volume of the star (r ^ 1). The density perturbation Sp appearing 
in equation (5) is calculated by taking the curl of both sides of equation (3) and matching the (()-components: 


dSp 

89 


r dr 
poRt, d"!* 


{V X [(V X B X 




Equations (5) and (6) are then solved to obtain e. 

We require the field to obey the following conditions: 


(i) the field is symmetric about the a-axis; 

(ii) the external field is current-free and is purely poloidal; 

(iii) the poloidal component of the field is continuous everywhere; 

(iv) the toroidal component of the field is confined to a toroidal region inside the star around the neutral curve; 

(v) the current is finite and continuous everywhere inside the star and vanishes at the surface of the star. 


( 6 ) 


These conditions are to be fulfilled by judicious choices of a and /3. 

Formally, the function /3(a) is arbitrary. In previous work (Reisenegger 2009; Mastrano et al. 2011; Akgiin et al. 2013; 
Mastrano, Lasky, & Melatos 2013), we take 


/3(a) 


(]a] - acfi for ]a] > ac, 
0 for ]a] < ac, 


( 7 ) 


where ac is the value of a that defines the last poloidal field line that closes inside the star. This relatively simple form of /3 
ensures that the toroidal field is confined to the equatorial torus bounded by the last closed poloidal field line. However, if 
the field line a(r, 9) — ac is not symmetric around some 9, this simple polynomial form for /3 leads to a discontinuity in Sp at 
the a = ac boundary (Mastrano, Lasky, & Melatos 2013). A worked example showing this explicitly for a quadrupolar field 
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is presented in the Appendix. For this reason, previous papers failed to construct multipolar fields with poloidal and toroidal 
components analytically (Mastrano et al. 2011; Mastrano, Lasky, & Melatos 2013). 


2.2 Non-barotropicity 

As noted above, we consider the star to be non-barotropic in this paper, i.e., we do not require a one-to-one relation between 
pressure and density, p — p{p)- Neutron star matter consists of multiple species (at least neutrons, protons, and electrons) 
which reach a stably-stratified, hydromagnetic equilibrium within a few Alfven time-scales (Pethick 1992; Reisenegger & 
Goldreich 1992; Reisenegger 2001; Riles 2013). This system is not in full chemical equilibrium: the relative abundances of the 
constituent particles change by weak nuclear interactions and diffusive processes. These processes have much longer time-scales 
than the Alfven time-scale (Hoyos, Reisenegger, & Valdivia 2008), with 

T 

-diff=1.7xl0«(^) yr, (9) 

respectively, where T is the temperature of the star. Between the Alfven and the weak nuclear time-scales, the star is in a 
hydromagnetic equilibrium state, in which the composition is not determined solely by the density or pressure, and density 
and pressure do not correspond one-to-one (Mastrano et al. 2011). Since most magnetars are ~ 1-10 kyr old (as inferred from 
spin and spin down), the calculations presented in this paper are readily applicable to them. In addition, note the strong 
inverse dependence on T; neutron stars cool down rapidly over ~ 10^ kyr (Yakovlev, Levenfish, & Shibanov 1999; Petrovich 
& Reisenegger 2010, 2011; Gonzalez-Jimenez, Petrovich, & Reisenegger 2014), so that, in practice, chemical equilibrium may 
never be reached. 


Tweak = 4.3 X 10“ 


2.3 Field-aligned coordinates 

In this section, we describe a method of generating a that does not lead to discontinuities in 5p while keeping equation (7) for 
/3(a). There are no cross terms between the poloidal and toroidal field components in the Lorentz force, so we can calculate 
the density perturbations 5pt and 5pp caused by the toroidal and poloidal field components separately and add them to arrive 
at the total density perturbation 5p. Substituting equation (1) into equation ( 6 ), we find 


and 


d$d(5pt d / ,da\ 

Bl dr de ^ de\ dr) 


dr 



( 10 ) 


R2 

^0 


d<l> d&pp 

d 

f l da 

1 d^a 


( 1 (3a \ 

dr do 

^ de ' 

( sin e dr 

sin 0 dr"^ 

r^de' 

/^sin 0 de J 


dr 


1 da 
sin 6 dO 


1 da 


sin 8 dr^ r^ d6 V sin 9 dd 


1 d^a 1 d 

+ 


( 11 ) 


for the two components, with F = (rsin0)“^ and P' = d/3/da. 

There are many valid ways to choose /3(a). However, numerical simulations (Braithwaite & Nordlund 2006; Braithwaite 
2009) favour configurations where the toroidal field is confined fully inside the star in a closed region around the neutral 
curve. We therefore seek a function a that, when substituted into /3(a), ensures continuity of 5pt at a = ac. (Note that 
5pp is continuous everywhere, since a is continuous everywhere). This task is made more tractable by first performing a 
coordinate transformation from spherical polar coordinates to those defined by the stream function a and some coordinate 
7 (which indicates angular position on the meridional field line defined by the coordinate a), i.e., {r,9,(j>) i—>■ (a, 7 , 0 ). The 
coordinates a and 7 must be chosen such that the Jacobian of this transformation is nonzero in the regions of interest and 
the transformation is invertible (do Carmo 2011). Furthermore, if 7 is chosen to satisfy djjdr = 0 everywhere, equation (10) 
can be rewritten as 


da dSptia,^) dj dSpAa,j) _ dF{a,j) ^^^^^,^^^dadj 

de da de d-i d-f ^ ’ dr de’ ^ ' 

with Jpj = poRtSptd^/dr. If a and 7 are chosen such that dajde, d^jde, and the entire right-hand side of equation (12) 
are all functions of a and 7 only, then equation ( 12 ) is a well-defined, first-order, linear, inhomogeneous partial differential 
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equation that permits a unique solution in closed form. We illustrate how to perform this construction through a specific 
example in Sec. 3. 

The mapping (r, 9, (j)) (a, 7 , (j>) is not unique. It depends not only on the particular magnetic field being investigated 

but also on one’s choice of 7 and is often cumbersome to write down. While the choice of a is obviously dictated by the poloidal 
magnetic field, the actual form of a can be complicated when one deals with a linear combination of many multipoles. For 
the relatively simple case of a dipole or dipole-plus-quadrupole field, we can choose (as we do in Sec. 3) a simple form for 7 , 
namely 7 = cos6, and equation (12) simplifies. However, when g{6) [the angular part of a(r,9)] is a combination of many 
powers of sin 9 and cos 9, it becomes difficult to choose 7 such that equation (12) is both well-defined and simple. 


3 WORKED EXAMPLE: DIPOLE-QUADRUPOLE, POLOIDAL-TOROIDAL FIELD 

In this section, we calculate Sp due to a mixed dipole and quadrupole magnetic field, with a toroidal component confined to 
the region bounded by the last poloidal field line that closes inside the star. The external field, to which the internal field is 
matched at r = 1, is given by (Mastrano, Lasky, & Melatos 2013) 


Rext — Rq 


2 cos 0 ^ fv(3cos^0— 1 ) 


©7- 



2k sin 9 cos 9 \ 
^4 J 



(13) 


where Bo sets the strength of the field overall and k is the dimensionless weight of the quadrupole component. 

The deformation caused by a field without an internal toroidal component matched to equation (13) has been calculated 
previously (Mastrano, Lasky, & Melatos 2013). Therefore, in this section, we focus on the deformation caused by the toroidal 
component. We begin by discussing the conditions that a must obey, then we derive a suitable (but not unique) form of a, 
i.e., one which generates an internal poloidal-toroidal magnetic field that matches Bext at r = 1, for which the Lorentz force 
is continuous everywhere inside the star. Lastly, we calculate Spt and hence Sp overall. 

Throughout the rest of this paper, to simplify analysis and to facilitate comparison with previous work, we use the 
gravitational potential, density, and pressure profiles used by Mastrano et al. (2011), namely 


P = Pc{l-r^), 

(14) 

P = P.(i-^ + 2^‘-y). 

(15) 


(16) 


where pc = 15 M*/( 87 rR*) and pc — 15GM*/(Ifirri?*) are the density and pressure at the origin, respectively. This simple, 
‘parabolic’ density profile has been shown to result in values of e that are ~ 5% of those found using the more realistic 
polytropic equation of state (Mastrano et al. 2011). 


3.1 Field-aligned coordinates 

The main problem with a magnetic field that is not symmetric around some 9 is making sure that Spt is continuous across 
a = Oc (Mastrano, Lasky, & Melatos 2013). An example showing why this does not occur automatically is worked out in the 
Appendix. As explained in Sec. 2.2, the issue can be circumvented by working in the field-aligned coordinates ( 0 , 7 , <()), where 
it is easier to choose a judiciously by inspection. 

Consider a stream function a that is a linear combination of dipole and quadrupole components, viz. 


a{r, 9) = /i(r) sin^ 9 + Kf 2 {r) sin^ 9cos9. (17) 

Field and current continuity at the surface of the star require that the radially dependent factors satisfy (Mastrano et al. 
2011; Mastrano, Lasky, & Melatos 2013) 


/i(l) = l = /2(l), (18) 

2/((l) = -2 = /^(l), (19) 

/('(I) - 2/1 (1) = 0 = fi ' il ) - 6/2(1), (20) 

where the prime indicates differentiation with respect to r. The first and second conditions ensure the continuity of Br and Bg 
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respectively. The third condition ensures current continuity (i.e., the current vanishes at the surface). Implicitly, the current 
is also well-behaved at the origin. 

The last closed internal poloidal field line is dehned by a = Oc, where Oc is the maximum of /i (1) sin^ 0 -|-k/ 2(1) sin^ 6 cos 9 
for 0^0^ 7r/2. In the case of a pure dipole, there Is only one neutral curve at 0 = 7r/2. In the case of a pure quadrupole, there 
are two neutral curves at 9 — cos“^(±-\/l/3). Mixed dipole-quadrupole cases may have one or two neutral curves, depending 
on the relative weights of the multipoles; see Fig. 2 of this paper and Fig. 3 of Mastrano, Lasky, & Melatos (2013) for some 
examples. For consistency and simplicity, we assume in this paper that the toroidal component only exists around the upper 
hemisphere’s neutral curve. 

Now, suppose we take 7 = 0036^. This choice is not unique, but it has the advantages that it is single-valued in the domain 
0 ^ ^ TT and that, trivially, we have d^/dr = 0. Next, we need to find suitable radial functions /i(r) and /2(r) that satisfy 

the above conditions (18)-(20), while bearing in mind that equation (12) is a well-dehned linear partial differential equation 
in a and 7, i.e., there are no ‘stray’ factors of r, sin^, or cos9 left. We choose /i(r) and /2(r) to be polynomials in r, which 
are simple and analytically tractable. We need at least three terms in each polynomial, so as to satisfy the three boundary 
conditions above. The implicit condition on the current at the origin is easily satished by choosing the starting powers of r 
high enough. 

In fact, even after applying these constraints, there are still an inhnite number of polynomials that meet our needs. We 
make the final choice through trial and error, to ensure that r(a, 7) is real, with 0 7^ r(a, 7) ^ 1 in the (0,7) domain of 
interest, and that it is expressible in closed form (for ease of calculation). We note that air" -|- a2r^" -I- is a cubic in r" 

and is guaranteed to have at least one real root. We note also that the lowest power of r for which the quadrupole’s current 
vanishes at the origin is n = 3, but /i(r) cannot have an r® term, if the dipole’s current is to vanish at the origin, so one 
needs n ^ 4. In addition, we add an arbitrary, extra term to the /i(r) polynomial, whose coefficient a is not hxed by any of 
the boundary conditions; it is left as a free parameter. Adjusting a changes the volume of the toroidal held region, a quantity 
of key physical importance, in a convenient fashion (Mastrano, Lasky, & Melatos 2013; Dall’Osso et al. 2014). Having chosen 
a value of a, one can then check for the existence of real roots in the domain 0 ^ r ^ 1; the existence is guaranteed for a — 0 
and for the values of a discussed in this paper. Now, choosing the lowest power of r in both /i(r) and /2(r) to be we 
Hnally arrive at 


= (W (S ■ + (i ■ 

/2(r) = I (35/ - 42r® -f 15r^^) . (22) 

8 

The choices leading to equations (21)-(22) are not unique, but they do guarantee a well-dehned inverse coordinate transform 
relation r(a,7) for all k by inverting a = /i(r)(l — 7^) -|- Kf2{r)~/{1 — 7^). We plot two examples of a(r,9) and 7(0) on the 
r-9 plane in Fig. 1 for k = 0 (i.e., a pure dipole, left-hand panel) and k = 2 (right-hand panel) with a = —2 hxed. Contours 
of a follow the poloidal held lines. The pure dipole case is symmetric around the equator, but the mixed case is not. As the 
right-hand panel shows, a mixed dipole-quadrupole case can have two neutral curves, one in each hemisphere. 

Eight example held-line plots with constant a and varying k (top four panels), and constant n and varying a (bottom 
four panels) are shown in Fig. 2. The top leftmost panel shows the k = 0 conhguration (i.e., a pure dipole). The other top 
panels show conhgurations with k = 0.2, 0.6, and 1.1. The three k ^ 0 plots show that a quadrupole component breaks the 
north-south symmetry of the held lines. Odd multipoles are north-south antisymmetric, while even multipoles are north-south 
symmetric; when odd and even multipoles are added together, the result is generally asymmetric (Mastrano, Lasky, & Melatos 
2013). The bottom panels show a conhguration with k. = 0.6 and different sizes of the toroidal region. As a increases, the 
toroidal held volume decreases. A pure dipole has one neutral curve at 9 = 7r/2, but the addition of a quadrupole shifts this 
original neutral curve northwards and introduces a new neutral curve into the southern hemisphere. As the bottom right 
panel of Fig. 2 shows, the new neutral curve also grows as the original toroidal region shrinks, even when k is kept constant. 
The general asymmetry of these composite helds, both internal and external, can potentially be linked to observables such as 
asymmetry in neutron star emission and burst signals (see Sec. 4). 

The helds in Fig. 2, generated by equations (17), (21), and (22), have smaller magnitudes at the origin than those analysed 
previously (Mastrano et al. 2011; Akgiin et al. 2013; Mastrano, Lasky, & Melatos 2013; Dall’Osso et al. 2014), as indicated by 
the sparseness of held lines around the origin. Indeed, some held lines seem to bend away from the origin, as in Fig. 2. This is 
a consequence of the higher powers of r in /i (r) and /2 (r) used in this paper compared to previous works. The detailed held 
behaviour at r = 0 is not a major concern here, as our primary focus is on stellar ellipticity, and the region near the origin 
contributes little to the moment of inertia. 
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Figure 1. Contour plots of field-aligned coordinates a{r,9) (solid curves) and 7 (r, 0) (dashed curves) on the r—9 plane [from equations 
(17), (21), and (22)] for a pure dipole (k = 0; left panel) and a mixed dipole-quadrupole (k = 2; right panel) with a = —2 fixed. 


3.2 Deformation 

The density perturbations due to the toroidal and poloidal magnetic components, Spt and 5pp, are readily calculated from 
equations (10) and (11) respectively. The total density perturbation, 5pp + Spt, is then substituted into equations (4) and (5) 
to calculate e. 

In Fig. 3, we plot e as a function of the parameter A (the ratio of internal poloidal field energy to total internal field 
energy; 0 ^ A ^ 1) for a canonical magnetar (M = 1.4 Mq, i?* = 10^ m. Bo = 5 x 10^° T) for fixed a and four different values 
of K. In Fig. 4, we plot e versus A for fixed k and four different values of a. In both figures, we also plot e(A) for a purely 
dipolar field with the same M, i?*, and Bo as the thick dashed curve (Mastrano et al. 2011). The plots show that prolateness 
increases as A decreases (i.e., as the toroidal magnetic field energy increases), similar to the results of Mastrano et al. (2011) 
and DairOsso et al. (2014). Figure 3 shows that, for a given ct, the curve shifts rightward and downward as the quadrupole 
field component strengthens, i.e., the star tends to be more prolate for the same A as k increases. In contrast, a star with a 
purely poloidal magnetic field becomes more oblate, as the quadrupole component strengthens (Mastrano, Lasky, & Melatos 
2013). The curves intersect at some points. For example, we find e{a = —5, k = 0) = e{a = —5,k = 0.2) at A « 0.6. If one 
were to detect a magnetar near these particular values of e and A, the analysis in this paper cannot distinguish between the 
two models. 

The curves in Fig. 4 show the same general trend: prolatenes increases as A decreases. The curves for ct = —10 (dashed) 
and CT = —15 (solid) indicate that the star becomes more prolate, as the region occupied by the toroidal field grows. However, 
the curves for ct = 10 (dashed-dotted) and cr = 15 (dotted) indicate that the star becomes more prolate, as the toroidal 
field region shrinks below some volume (for this particular example, the star becomes more prolate again for a > 2). This 
counterintuitive trend occurs, because a second neutral curve (and the corresponding region with low poloidal field strength) 
emerges inside the star in configurations with high and positive a. The weak poloidal field around the second neutral curve 
makes it easier to deform the star into a prolate shape. The curves intersect at some points, like in Fig. 3. For example, we 
find e{a = —15,k = 0.2) = e{a = —10,«; = 0.2) at A = 0.58, e{a = 10, k = 0.2) = e{a — 15, k = 0.2) at A = 0.78, and 
e{a = —10, K = 0.2) = e{a = 15, k = 0.2) at A = 0.86. 

Yoshida (2013) calculated e for a pure dipole without taking the Cowling approximation. His values of |e| are ~ 2 times 
those found by Mastrano et al. (2011). Furthermore, Yoshida (2013) found that, if the zeroth-order density profile p(r) features 
an off-centred maximum (a somewhat artificial situation), taking the Cowling approximation results in the opposite shape 
from that predicted by the full perturbation analysis (i.e., a prolate star where the full perturbation analysis predicts an 
oblate star and vice versa). A thorough, full-perturbation calculation without the Cowling approximation is beyond the scope 
of this paper but is certainly worth undertaking in future. 
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Figure 2. Magnetic field lines of a mixed dipole-quadrupole configuration, as a function of toroidal volume (related to a) and quadrupole 
coefficient k. The internal field is described by equations (1), (2), and (17)—(22). The external field is described by equation (13). Each 
panel is labelled by (ct, k). Top panels: a = 0.5; k = 0, 0.2, 0.6, and 1.1 (left to right). Bottom panels: k = 0.6; a = —20, —10, 0, and 
4 (left to right). In each panel, the shading defines the volume to which the toroidal magnetic field component is confined. The dashed 
semicircle represents the stellar surface. 


3.3 Magnetic moment of inertia 

In Sec. 3.2, we neglect the direct contribution to e from the magnetic energy density in the component of the stress-energy 
tensor.^ Most recent works on this subject [e.g., Haskell et al. (2008); Dall’Osso et al. (2014)] neglect this contribution as well. 
It is given by (Thorne 1980) 


es = [ [ dr dO sm6{l — 3cos^ 9) ^ , (23) 

lo Jo Jo 2^0 

which is essentially the same as equation (13) of Mastrano et al. (2011), with Sp replaced by l2po<?■ The external field 
contribution to es is comparable to equation (23), because externally we have Bext oc [equation (13)] and dr r'^B^ is 
dominated by the lower terminal. For the dipole-plus-quadrupole example in Sec. 3.2 with a = —5 and k = 1.1 (corresponding 
to the top rightmost panel of Fig. 2 and the thin dotted curve in Fig. 3), we find 


€b = -1.4 X 10” 


Bo 


5 X IQio T 


M* 


1.4 Me 


10^ m J 



0.847\ 
A )' 


(24) 


^ We thank the reviewer for bringing this issue to our attention. 
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Figure 3. Mass ellipticity e (in units of 10“®) versus A (the ratio of internal poloidal magnetic field energy to total internal magnetic 
field energy), for a = —5 fixed and «: = 0 (thin solid curve), 0.2 (thin dashed curve), 0.6 (thin dashed-dotted curve), 1.1 (thin dotted 
curve), for a canonical magnetar {M = 1.4 Mq, i?* = 10“^ m, Bq = 5 x 10^® T). The result for a pure dipole (Mastrano et al. 2011) is 
included for comparison (thick dashed curve). 



A 

Figure 4. Mass ellipticity e (in units of 10“®) versus A (the ratio of internal poloidal magnetic field energy to total internal magnetic 
field energy), for k, = 0.2 fixed and a = —15 (thin solid curve), —10 (thin dashed curve), 10 (thin dashed-dotted curve), 15 (thin dotted 
curve, for a canonical magnetar (M = 1.4 Mq, i?* = lO'^ m, Bq = 5 x 10^® T). The result for a pure dipole (Mastrano et al. 2011) is 
included for comparison (thick dashed curve). 
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A 


Figure 5. Mass ellipticity e (in units of 10“®) versus A (the ratio of internal poloidal magnetic field energy to total internal magnetic 
field energy), for a dipole-plus-quadrupole field configuration, with a = —5 and k = 1.1, where we take into account (solid curve) and 
neglect (dashed curve) the moment-of-inertia tensor contribution of the magnetic field, for a canonical magnetar (M = 1.4 Mq, R* = 10"^ 
m, Bq = 5 X 10^® T). The two curves nearly overlap. 


Note that es scales differently with R, and M* [compared to the 5p contribution to e, e.g., given by equation (25) of this 
paper or equation (17) of Mastrano et al. (2011)], because es is not coupled to V4> [given by equation (16) in this paper for 
the parabolic density profile]. 

We plot e versus A for a canonical magnetar (M = 1.4 Mq, Rt = 10 “* m. So = 5 x 10 ^° T), with a = —5 and k = 1.1 
in Fig. 5. The figure shows two nearly overlapping curves. The dashed (solid) curve corresponds to excluding (including) the 
magnetic field contribution. Our calculation confirms that, at worst, the effect shifts e by a factor of ~ 2%. This is because 
one has 5pj{B^ jpoc^) ~ c^(dr/d$) ^ 1. 

The focus of this methods paper is on the star-deforming effect of the magnetic field. The uncertainty introduced by 
neglecting es is smaller than that potentially introduced by the commonly-taken Cowling approximation (Yoshida 2013), nor 
does it change the general qualitative behaviour of the results. 


4 ASTROPHYSICAL EXAMPLES 

While this paper is mainly a ‘methods’ paper, we discuss briefly, by way of illustration, two examples of how one can apply 
the method to model astrophysical phenomena in the magnetars SGR 0418-1-5729 (Giiver, Ozel F., & G6gu§ 2008; Giiver, 
G6gu§, & Ozel 2011; Tiengo et al. 2013) and 4U 0142+61 (Makishima et al. 2014). 

4.1 SGR 0418+5729 

SGR 0418+5729 has an inferred dipole field strength of < 7.6 x 10® T (Rea et al. 2010), but an analysis of the X-ray spectrum 
by Giiver, G6gu§, & Ozel (2011) concluded that a surface field strength of 10^° T fits the data best. Giiver, G6gu§, & Ozel 
(2011) suggested that the surface field contains higher-order multipole(s) (which fall away with r faster than the dipole) to 
account for the discrepancy. If we take Bq = 3.8 x 10® T, corresponding to a dipole field strength of 7.6 x 10® T at the polar 
surface, and take k = 13, corresponding to total (dipole plus quadrupole) field strength of 10^® T at the polar surface, then 
we find 


e = 1.7 X 10~ 


Bo 


5 X lOio T 


M* 


R* 


A - 1 + 7 X 10” 


^1.4Moy VlOkm; V A J’ 

for CT = 0. We choose to write e in terms of Bq (which can, in principle, be inferred from spin and spin down) rather than 
(say) the internal magnetic field strength (maximum or volume-averaged) in order to relate e to an observable quantity, which 
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Table 1. Sample models of 4U 0142+61 for t = —1.6 X 10“^, the elllpticity required by interpreting X-ray pulse 
modulation as precession (Makishima et al. 2014), for some different values of k, with cr = —0.1 fixed. We show 
K (first column), total polar surface field strength Bpoie (second column), the required poloidal-to-total field 
energy ratio (third column), and the maximum toroidal field strength Bt (fourth column). 


K 

Bpoie (10i2 T) 

A, 

Bt (10i2 T) 

0.78 

4.7 X 10-2 

2.8 X 10-1 

16 

1 

5.2 X 10-2 

9.2 X 10-1 

9.5 

2 

7.8 X 10-2 

2.7 X 10-3 

7.1 

5 

1.6 X 10-1 

6.2 X 10-3 

9.3 

10 

2.9 X 10-1 

3.8 X 10-2 

6.9 


does not itself depend on A. Furthermore, this facilitates comparison with other works (Mastrano et al. 2011; Dall’Osso et al. 
2014). 

The above model does not need much toroidal component to deform into a prolate shape; we obtain e ^ 0 for 1 — A > 
7 X 10“^. If SGR 0418+5729 has no toroidal field component, we obtain e = 6.8 x lO“®(M*/1.4M0)“^(i7«/lO km)^, which is 
too small to generate gravitational waves detectable by current-generation interferometers (Mastrano et al. 2011; Riles 2013). 
However, A = 0.9 is enough to result in e = — 1O“®(M«/1.4M0)“^(R*/1O km)^, in contrast with a purely dipolar neutron 
star, which only becomes significantly prolate for A ^ 0.4 (Mastrano et al. 2011). Thus, in objects like SGR 0418+5729, e 
depends strongly on field configuration, characterized by A, as well as field magnitude. Future gravitational-wave upper limits 
from next-generation detectors will constrain A better in the context of this model. Because of the low spin frequencies of 
magnetars, e.g., 0.11 Hz for SGR 0418+5729 (G6gu§, Woods & Kouveliotou 2009), the best limits will come from the Einstein 
Telescope, where Newtonian noise suppression will extend observations down to the sub-Hertz regime (Bennett, van Eysden, 
& Melatos 2010). Present-day limits on e complement historical limits based on fast-spinning birth scenarios (Thompson, 
Ghang, & Quataert 2004; Dall’Osso & Stella 2007; Dall’Osso, Shore, & Stella 2009; Melatos & Priymak 2014). 


4.2 4U 0142+61 

The magnetar 4U 0142+61 has an inferred dipole field strength of 1.3 x 10^° T (Gavriil & Kaspi 2002; Giiver, Ozel F., & 
G6gu§ 2008). Phase modulations of its hard X-ray pulses, observed by the Suzaku satellite (Enoto et al. 2011; Makishima 
et al. 2014), have been interpreted by Makishima et al. (2014) as evidence of free precession, indicating a prolate star with 
|e| ~ 1.6 X 10“^ and, hence, a maximum internal toroidal field Bt ~ 10^^ T. 

As in the previous subsection, we assume that the inferred dipole field is not the whole story: the actual field contains 
a quadrupole element, which cannot be inferred from spin down, and an internal toroidal field, whose presence can only be 
detected indirectly via its effect on e. From equation (13), the polar surface field strengths of the dipole and quadrupole 
components are Bdip = 2Bo and Squad = 2kBo respectively. If we infer that Sdip = 1-3 x 10^° T from spin down, we can 
calculate the values of A required to obtain e = —1.6 x 10“^, as required by the precession interpretation, for different values 
of K (with a = —0.1 in all cases). We present the cases of «: = 0.78, 1, 2, 5, and 10 in Table 1. Each row shows k, Bpoie, which is 
the total polar surface field strength, A^, which is the value of A required to impart e = —1.6 x 10~^ to the star, and Bt, which 
is the maximum strength of the toroidal field. Note that k = 0.78 is chosen because this corresponds to Bpoie = 4.7 x 10^° T, 
the value obtained by Giiver, Ozel F., & G6gu§ (2008) when they analysed the X-ray pulses from 4U 0142+61. 

For the cases investigated, we find values of Bt that are ~ 10 times stronger than calculated by Makishima et al. (2014). 
This, however, can be counteracted by reducing a and, hence, enlarging the volume occupied by the toroidal field. Table 1 
shows that, as the quadrupole component strengthens, Ae increases, i.e., a more quadrupolar field needs less toroidal field 
energy to deform the star. However, it is harder to say with confidence what trend Bt follows, e.g., Bt decreases between 
K = 0.78 and 2, then increases between k — 2 and 5, then decreases again between «: = 5 and 10. This is because k affects the 
toroidal field’s volume as well, although to a lesser degree than a. Thus, for example, we find Ae(K = 2) < Ae(«: = 5), i.e., the 
K = 5 case requires less toroidal field energy to deform the star to the desired e. However, we also find Bt{K = 2) < Bt{K = 5), 
because the toroidal field’s volume for k = 5 is smaller than that for k = 2. 

Note that we do not claim the magnetar 4U 0142+61 to be in any of the magnetic configurations explored in this 
subsection. We present these examples simply to demonstrate what our analytical approach can accommodate and analyse. 


5 CONCLUSION 

In this short methods paper, we generalise previous calculations (Mastrano et al. 2011; Mastrano, Lasky, & Melatos 2013) 
to show how any multipolar stellar magnetic field containing both poloidal and toroidal axisymmetric components can be 
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constructed analytically to satisfy boundary conditions (zero surface currents and zero toroidal field outside the equatorial 
torus) motivated physically and by numerical simulations in the literature. We also present a worked example of a dipole- 
plus-quadrupole, poloidal-plus-toroidal field and calculate e versus A for some representative combinations of the parameters 
a (which controls the volume occupied by the toroidal field) and k (which controls the weight of the quadrupole component). 
The star tends to be more prolate as k increases (Fig. 3). In the fixed k = 0.2 example shown in Fig. 4, for a < 2, the star 
deforms into a more prolate shape as a decreases. For cr > 2, increasing a makes the star more prolate, as the internal poloidal 
held component weakens. For a given A, i.e., for a given toroidal held energy, smaller a means greater toroidal held strength. 

In Sec. 4, we briehy discuss two possible astrophysical applications of our analysis to magnetars. SGR 0418+5729 is 
interesting because there is a mismatch between the surface held strengths inferred from spin down (Rea et al. 2010) and from 
X-ray analysis (Giiver, G6gu§, & Ozel 2011). We show how, if the star’s magnetic held is mostly quadrupolar [as suggested by 
Giiver, G6gu§, & Ozel (2011)], the star is likely to be prolate. For k = 13, we hnd 1 — A ^ 7 x 10“^ for e ^ 0 (i.e., the toroidal 
held only needs to contribute at least 0.07% to the total held energy to deform the star into a prolate shape). Present-day 
and historical (at birth) upper limits on gravitational wave emissions can thus be used to infer the internal held conhguration 
and toroidal held strength. 4U 0142+61 is interesting, because Suzaku observations suggest that it undergoes free precession, 
from which an upper limit for e can be calculated (Makishima et al. 2014). We present several possible conhgurations of this 
magnetar’s held and calculate the lower limit on A implied by the upper limit on e. We hnd that, to obtain jej ~ 10“^, we 
need toroidal helds with maximum strength Bt ~ 10^® T, ten times greater than inferred by Makishima et al. (2014). The 
two applications are by no means exhaustive; they simply illustrate the potential of the method. 

Throughout this paper, we do not discuss the stability of these dipole-plus-quadrupole, poloidal-plus-toroidal held conhg¬ 
urations. Such calculations are reserved for future work. The results of Akgun et al. (2013) are not directly applicable to the 
composite multipolar helds discussed here. To repeat their analysis, one must hrst identify the regions of greatest instability in 
the star. This is more difficult than for the pure dipole held, because the held is no longer north-south symmetric, because the 
toroidal held is ohset from the equator, and because there are new regions of low held strength inside the star (that is, other 
than the neutral curve). It may be easier to test the stability of these composite helds numerically rather than analytically, 
by evolving the held in a time-dependent magnetohydrodynamic simulation. We defer this calculation to a future paper. 
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APPENDIX A: AN EXAMPLE OE A DISCONTINUOUS Spt ARISING FROM A POOR CHOICE OF 

fir) 

In this short Appendix, we present an example of what happens when one uses the simple form of Pia) given in equation (7) 
without choosing the polynomial /(r) judiciously. We start with a pure quadrupole field. 


ct = f 2 (r) sin^ 0 cos 6, 


(Al) 


/3(a) 


(]a] - ac)^ for ]a] > a^, 
0 for ]a] < ac. 


(A2) 


with ac = 2\/3/9. In this section, we choose the radial function /2(r) = 21(r® — |r^ + |r®), as chosen by Mastrano, Lasky, 
& Melatos (2013). This /(r) is a natural choice, because it fulfills the boundary conditions (i)-(v) listed in Sec. 2.1, because 
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Figure Al. Density perturbation due to the toroidal field component, 5pt, versus colatitude 9 for r = 0.7, for a quadrupole where the 
field is described by equations (A1)-(A2) (left-hand panel) and for a pure dipole (right-hand panel). 


it has the minimum number of terms, and because its first power, r®, is the lowest power of r that guarantees the current 
vanishes at the origin. However, this choice of /(r) does not allow us to perform an invertible coordinate transformation 
(r, 0,(^) (a,7, (/>), and therefore does not permit a unique solution in closed form for equation (12). 

Substitution of equations (A1)-(A2) into equation (10) leads to 


UoR* 


R2 


d<F dSpt 

0 

2{a — Oc)^ Oa 

0 

2{a — Oc)^ Oa 

dr 09 

^ 09 

sin^ 9 Or 

Or 

r^ sin^ 9 09 


(A3) 


Integrating equation (A3) with respect to 0 should give Spt, plus some integrating constant g{r), to be calculated by matching 
Spt = 0 at a = Oc. We plot 5pt versus 9 for r — 0.7 for the quadrupole in the left-hand panel of Fig. Al as an example. We obtain 
Spt{r — 0.7, 9 — 9i) ^ Spt{r = 0.7, 9 — @2), with 9i and 92 being the coordinates where a{r — 0.7, 9) = a^. For comparison, 
the right-hand panel of Fig. Al shows 5pt{r = 0.7) for a pure dipole, showing Spt{r — 0.7,9 = 9i) = Spt{r = 0.7, 6^ = @2) 
(Mastrano et al. 2011). It is therefore impossible to ensure the continuity of Spt at a = for a quadrupole (unlike for the 
dipole) using only an arbitrary function of r. This problem vanishes when the method described in Sections 2.2 and 3 is 
applied. 
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